In this document I compare gene-snp pairs with strong heterogeneity:
knitr::opts_chunk$set(cache=TRUE)
covmat=readRDS("../withoutbrain/covmatnobrain.rds")
pis=readRDS("../withoutbrain/pisnobrain.rds")$pihat
z.stat=read.table("../Data/maxz.txt")
b.stat=read.table("../Data/maxbetahats.txt")
pi.mat=matrix(pis,ncol=44,byrow=T)
lfsr.mash=read.table("../withoutbrain/nobrainlfsr.txt")[,-1]
pm.mash=read.table("../withoutbrain//nobrainposterior.means.txt")[,-1]
pm.mash.brain=read.table("../Data/Aug13withEDposterior.means.txt")[,-1]
standard.error=read.table("../Data/standard.error.txt")
pm.mash.beta.brain=pm.mash.brain*standard.error
pm.mash.beta=pm.mash*standard.error[,-c(7:16)]
pm.ash=read.table("../Comparison with Univariate Ash/univariate.ash.pm.txt")[,-c(7:16)]
lfsr.ash=read.table("../Comparison with Univariate Ash/univariate.ash.lfsr.txt")[,-c(7:16)]
pm.ash.beta=pm.ash*standard.error[,-c(7:16)]
colnames(lfsr.mash)=colnames(pm.mash)=colnames(pm.mash.beta)=colnames(z.stat)[-c(7:16)]
het.norm=function(effectsize){
t(apply(effectsize,1,function(x){
x/x[which.max(abs(x))]
}))}
het.func=function(normdat,threshold){
apply(abs(normdat),1,function(x){sum(x>threshold)})}
matrix.ash.norm=het.norm(effectsize = pm.mash.beta)
uni.ash.norm=het.norm(effectsize = pm.ash.beta)
matrix.ash.norm.brain=het.norm(effectsize = pm.mash.beta.brain)
het.genes.no.brain=het.func(matrix.ash.norm,0.5)
het.genes.with.brain=het.func(matrix.ash.norm.brain,0.5)
Let’s looks at how heterogeneity behaves:
library(qtlcharts)
iplotCurves(curveMatrix = pm.mash.brain[het.genes.with.brain>40,])
## Set screen size to height=700 x width=1000
iplotCurves(curveMatrix = z.stat[het.genes.with.brain>40,])
Now, we look at singletons:
iplotCurves(curveMatrix = pm.mash.brain[het.genes.with.brain<2,],chartOpts=list(curves_xlab="Tissue", curves_ylab="PMMASHBrain"))
iplotCurves(curveMatrix = z.stat[het.genes.with.brain<2,],chartOpts=list(curves_xlab="Tissue", curves_ylab="Z.statBrain"))
Let’s consider what happens when we remove brains:
iplotCurves(curveMatrix = pm.mash.beta[het.genes.no.brain >40,])
iplotCurves(curveMatrix = b.stat[het.genes.no.brain>40,-c(7:16)])
Now, we look at singletons:
iplotCurves(curveMatrix = pm.mash.beta[het.genes.no.brain <2,])
iplotCurves(curveMatrix = b.stat[het.genes.no.brain<2,-c(7:16)])